function [parties,voters]=Politics_Laver_2D(n_parties,name,name_abbr,phi_1,phi_2,phi_3, ...
                                          gamma,kappa,delta_val,delta_val_change,val_start,val_1,val_change,val_2,n_it,t_start,pos_init, ...
                                          n_voters,Density_CF,Density_PL,CF_vec,PL_vec,space_size_1,space_size_2,correlation, ...
                                          max_it,polls,show_num,show_graph,create_mov,color_vec,n)
   %% Declare and initialize parties (candidates) and voters structures/parameters
   
   parties.phi_1=zeros(n_parties,max_it);           %Parties strategy
   parties.phi_2=zeros(n_parties,max_it);
   parties.phi_3=zeros(n_parties,max_it);
   for i=1:n_parties
      parties.phi_1(i,1:max_it)=phi_1(i);
      parties.phi_2(i,1:max_it)=phi_2(i);
      parties.phi_3(i,1:max_it)=phi_3(i);
   end
   parties.name=name;                               %Name to refer to in the figures
   parties.number=1:n_parties;                      %Parties ID
   parties.psize=zeros(n_parties,max_it);           %Parties vote share
   parties.xcor=zeros(n_parties,max_it);            %Parties positions in 2D space
   parties.ycor=zeros(n_parties,max_it);
   parties.heading=zeros(n_parties,1);              %Parties heading direction
   c3=[1.0 1.0 1.0];                                %Color in figure
   c2=[0.0 0.3 1.0];
   c1=[1.0 0.0 0.0];
   parties.color(:,:,1)=parties.phi_3*c1(1)+parties.phi_2*c2(1)+parties.phi_1*c3(1);
   parties.color(:,:,2)=parties.phi_3*c1(2)+parties.phi_2*c2(2)+parties.phi_1*c3(2);
   parties.color(:,:,3)=parties.phi_3*c1(3)+parties.phi_2*c2(3)+parties.phi_1*c3(3);
   parties.shape='^^^^^^^^';                        %Shape in figure
   parties.Rp=zeros(max_it,1);                      %Representativeness of whole political field
   parties.valence=repmat(val_1',1,max_it).*zeros(n_parties,max_it); %Parties valence
   for i=1:n_parties
      parties.valence(i,t_start(i):t_start(i)+delta_val-1)=linspace(val_start(i),val_1(i),delta_val);
      parties.valence(i,t_start(i)+delta_val:end)=val_1(i);
   end
   
   %Valence change every delta_val_change weeks (starting from week 17)
   timeframe(1:n_it+1)=17:delta_val_change:max_it;
   for t=1:n_it
       if t>=2
          i=val_change(t-1);
          if i~=0
             j=timeframe(t);
             parties.valence(i,j:j+delta_val_change)=linspace(parties.valence(i,j),val_1(i),delta_val_change+1);  
             parties.valence(i,j+delta_val_change+1:end)=val_1(i);
          end
       end
       i=val_change(t);
       if i~=0
          j=timeframe(t);
          parties.valence(i,j:j+delta_val_change)=linspace(parties.valence(i,j),val_1(i)+val_2(i),delta_val_change+1);  
          parties.valence(i,j+delta_val_change+1:end)=val_1(i)+val_2(i);
       end
   end

   %Voters' structure
   voters.closest_party=zeros(n_voters,max_it);
   voters.distance=zeros(n_voters,max_it);
   voters.Rp=zeros(max_it,1);
   voters.xcor=zeros(n_voters,1);
   voters.ycor=zeros(n_voters,1);
   
   %Set up parties and voters distributions and intiial positions
   [parties,voters,voters_distribution]=...
                          setup(parties,voters,n_parties,n_voters,Density_CF,Density_PL,CF_vec,PL_vec,...
                          pos_init,correlation);
        
   %Finds if any candidate is out of the race (either not yet candidate or dropped out)
   %by looking if the polls take them into account.
   candidate_in=find(polls(1,:)~=0);
   
   %Initialize voters preferences
   [parties,voters,largest_party]=update(parties,voters,n_parties,n_voters,1,candidate_in);

   %Show initialization results
   if strcmp(show_graph,'yes')
      h1=figure(1);
      h2=figure(2);
      show_results(parties,n_parties,n_voters,space_size_1,show_num,1,...
                   voters_distribution,CF_vec,PL_vec,candidate_in,create_mov,n,name_abbr,color_vec);
   end

   %------------------Time loop----------------
   for it=2:max_it

       %Update position of each party according to their strategy
       candidate_in=find(polls(it,:)~=0);
       for i=1:n_parties
           
           if sum(i==candidate_in)>=1

              %Move only if the candidate is below the satisfying threshold
              if parties.psize(i,it-1)/n_voters<kappa(i)
                 [parties]=hunt(parties,gamma(i),voters,i,pos_init(i,:),n_voters,it);
              else
                 parties.xcor(i,it)=parties.xcor(i,it-1);
                 parties.ycor(i,it)=parties.ycor(i,it-1); 
              end

           else
               
              parties.xcor(i,it)=parties.xcor(i,it-1);
              parties.ycor(i,it)=parties.ycor(i,it-1);
               
           end
           
       end   

       %Update voters preferences
       [parties,voters,largest_party]=update(parties,voters,n_parties,n_voters,it,candidate_in);

       %Show results
       if strcmp(show_graph,'yes')
          show_results(parties,n_parties,n_voters,space_size_1,show_num,it,...
                   voters_distribution,CF_vec,PL_vec,candidate_in,create_mov,n,name_abbr,color_vec);
       end

   end
   
   if strcmp(show_graph,'yes')
      close(h1)
      close(h2)
   end

   return  
     
end


%% FUNCTIONS

%-----------------------------SETUP FUNCTION-------------------------------
function [parties,voters,voters_distribution]=...
         setup(parties,voters,n_parties,n_voters,Density_CF,Density_PL,CF_vec,PL_vec,...
               pos_init,correlation)

   %Initialize position, strategy, heading and color of parties
   parties.xcor(:,1)=pos_init(:,1);
   parties.ycor(:,1)=pos_init(:,2);
   parties.heading=random('Uniform',0,360,[n_parties,1]);
   
   %Empirical distribution with correlation between two dimensions
   N_vec=size(CF_vec,2);
   rand_set=randpdf_2D(Density_CF,Density_PL,CF_vec,PL_vec,[n_voters,1]);
   voters.xcor=rand_set(:,1);
   voters.ycor=correlation*rand_set(:,1)+sqrt(1-correlation^2)*rand_set(:,2);
   voters_distribution=zeros(N_vec,N_vec);
   for i=1:n_voters
      ind_1 = ceil(voters.xcor(i)/2*N_vec);
      ind_2 = ceil(voters.ycor(i)/3*N_vec);
      if ind_1<1
         ind_1=1;
      end
      if ind_1>N_vec
         ind_1=N_vec;
      end
      if ind_2<1
         ind_2=1;
      end
      if ind_2>N_vec
         ind_2=N_vec;
      end
      voters_distribution(ind_1,ind_2) = voters_distribution(ind_1,ind_2) + 1;    
   end
      
   %disp(['Percentage of voters pro-choice: ' num2str(size(find(voters.ycor<=0.5*space_size_2),1)/n_voters*100,3) '%']);
   %disp(['Correlation between two dimensions is ' num2str(corr(voters.xcor,voters.ycor))]);
   
   return
   
end

%------------------------UPDATE FUNCTION-----------------------------
function [parties,voters,largest_party]=update(parties,voters,n_parties,n_voters,it,candidate_in)

   %Number of candidates to consider
   n_candidate_in=length(candidate_in);
   
   %Euclidean norm to get the distance squares
   norm = -( (repmat(parties.xcor(candidate_in,it),1,n_voters)-repmat(voters.xcor',n_candidate_in,1)).^2 + ...
             (repmat(parties.ycor(candidate_in,it),1,n_voters)-repmat(voters.ycor',n_candidate_in,1)).^2 )./...
              repmat(parties.valence(candidate_in,it),1,n_voters).^2;
 
   %Find closest party
   [voters.distance(:,it),voters.closest_party(:,it)] = max(norm,[],1);
   voters.closest_party(:,it) = candidate_in(voters.closest_party(:,it));
      
   %Compute representativeness of configuration
   voters.Rp(it) = -sum(voters.distance(:,it).^2)/n_voters;
   
   %Count number of supporters for each party
   for i=1:n_parties
       parties.psize(i,it)=length(find(voters.closest_party(:,it)==i));
   end
   
   %Find largest party
   [~,largest_party]=max(parties.psize(:,it));

   return

end

%------------------------HUNT FUNCTION-----------------------------
function [parties]=hunt(parties,gamma,voters,party_number,pos_init,n_voters,it)
 
   %Compute vote share, distance to ideal policy point and distance to
   %voter's centroid
   if it~=2

      vote_share_previous = parties.psize(party_number,it-2)/n_voters;
      vote_share_current = parties.psize(party_number,it-1)/n_voters;
      dist_policy_previous = sqrt((parties.xcor(party_number,it-2)-pos_init(1,1))^2 + (parties.ycor(party_number,it-2)-pos_init(1,2))^2);
      dist_policy_current = sqrt((parties.xcor(party_number,it-1)-pos_init(1,1))^2 + (parties.ycor(party_number,it-1)-pos_init(1,2))^2);
      supporters_current = find(voters.closest_party(:,it-1)==party_number);

      if isempty(supporters_current)==0
          mean_x_supp = mean(voters.xcor(supporters_current));
          mean_y_supp = mean(voters.ycor(supporters_current));
          dist_supporters_previous = sqrt( (parties.xcor(party_number,it-2) - mean_x_supp)^2  + ...
                                           (parties.ycor(party_number,it-2) - mean_y_supp)^2  );
          dist_supporters_current = sqrt( (parties.xcor(party_number,it-1) - mean_x_supp)^2  + ...
                                          (parties.ycor(party_number,it-1) - mean_y_supp)^2  );
      else
          dist_supporters_previous = 0;
          dist_supporters_current = 0;
      end

      %Compute utility function
      utility_previous = parties.phi_3(party_number,it-2)*vote_share_previous ...
                       - parties.phi_1(party_number,it-2)*dist_policy_previous^2 ...
                       - parties.phi_2(party_number,it-2)*dist_supporters_previous^2;
      utility_current = parties.phi_3(party_number,it-1)*vote_share_current ...
                      - parties.phi_1(party_number,it-1)*dist_policy_current^2 ...
                      - parties.phi_2(party_number,it-1)*dist_supporters_current^2;

      %If utility function has decreased (or kept the same), reverse direction            
      if utility_current <= utility_previous
         parties.heading(party_number)=parties.heading(party_number)-180+random('Uniform',-90,90);
      end

   end

   %Go to new point
   parties.xcor(party_number,it)=parties.xcor(party_number,it-1)+gamma*cos(pi/180*parties.heading(party_number));
   parties.ycor(party_number,it)=parties.ycor(party_number,it-1)+gamma*sin(pi/180*parties.heading(party_number));
   
   return

end

%------------------------SHOW_RESULTS FUNCTION-----------------------------
function show_results(parties,n_parties,n_voters,space_size,show_num,it,...
                      voters_distribution,CF_vec,PL_vec,candidate_in,create_mov,n,name_abbr,color_vec)


   %End simulation if figure has been closed
   f1 = get(groot);
   numfig=size(f1.Children);
   if numfig<2
      error(['Simulation ended at step ' num2str(it)]); 
   end
   
   %Plot position of voters and parties
   f=figure(1);
   if it==1
      clf;
      pcolor(CF_vec,PL_vec,voters_distribution');shading interp;hold on;
      set(gcf,'Color','w','Units','normalized','Position',[0.5 0 0.5 1]);
      set(gca,'Linewidth',2,'Box','on','FontSize',18,'TickLabelInterpreter',...
       'latex','XLim',[0 2],'YLim',[0 2]);
      xlabel('CF','Interpreter','latex');
      ylabel('PL','Interpreter','latex');
      axis square;
   end
   
   %Create orientated patch for each party
   h1=gobjects(n_parties,1); 
   h2=gobjects(n_parties,1); 
   h3=gobjects(n_parties,1); 
   s_patch=space_size/30;
   [x_patch,y_patch]=create_patch(parties.xcor(:,it),parties.ycor(:,it),parties.heading,...
                                  parties.shape,n_parties,s_patch);
         
   for i=candidate_in
      [h1(i)]=patch(x_patch(i,:),y_patch(i,:),squeeze(parties.color(i,it,1:3))');hold on; 
   end
   if it>=2
      for i=candidate_in
         line([parties.xcor(i,it-1) parties.xcor(i,it)],[parties.ycor(i,it-1) parties.ycor(i,it)],...
                'Color',squeeze(parties.color(i,it,1:3)),'Linewidth',1.0);hold on;
      end
   end
   if strcmp(show_num,'yes')
      for i=candidate_in
          [h2(i)]=text(parties.xcor(i,it)+0.8*s_patch,parties.ycor(i,it),...
               num2str(parties.psize(i,it)/n_voters*100,3),'Interpreter','latex',...
               'FontSize',12,'LineStyle','none','Color','w');hold on;
          [h3(i)]=text(parties.xcor(i,it)-1.5*s_patch,parties.ycor(i,it),...
               parties.name(i),'Interpreter','latex',...
               'FontSize',12,'LineStyle','none','Color','w');hold on;
      end
   end
   h=[h1,h2,h3];
   
   %Create movie
   if strcmp(create_mov,'yes')==1
      frame = getframe(f); 
      im = frame2im(frame); 
      [imind,cm] = rgb2ind(im,256); 
      % Write to the GIF File 
      if it == 1 
          imwrite(imind,cm,['Movie.gif'],'gif', 'Loopcount',inf,'DelayTime',0.2); 
      else 
          imwrite(imind,cm,['Movie.gif'],'gif','WriteMode','append','DelayTime',0.2); 
      end 
   end

   %Plot number of supporters per party
   figure(2)
   if it==1
      clf;
      set(gcf,'Color','w','Units','normalized','Position',[0 0 0.5 1]);
      set(gca,'Linewidth',2,'Box','on','FontSize',18,'TickLabelInterpreter',...
         'latex','XLim',[1 50]);
      xlabel('Iterations','Interpreter','latex');
      ylabel('Parties size [\%]','Interpreter','latex');    
   elseif it>=2
      for i=1:n_parties
         col=squeeze(color_vec(i,:));
         line([it-1 it],[parties.psize(i,it-1)/n_voters*100 parties.psize(i,it)/n_voters*100],...
             'Color',col,'Linewidth',2);hold on;
         scatter(it-1,parties.psize(i,it-1)/n_voters*100,'MarkerFaceColor','w',...
                 'MarkerEdgeColor',col,'Marker','o','SizeData',44);hold on;
      end 
   end
   
   drawnow
   %Delete objects
   delete(h);
   
   return

end

%------------------------CREATE_PATCH FUNCTION-----------------------------
function [x_patch,y_patch]=create_patch(xcor,ycor,heading,shape,n_parties,s_patch)

    %Create patch depending on party type and position
    x_patch=zeros(n_parties,4);
    y_patch=zeros(n_parties,4);
    for i=1:n_parties
       angle_rot=0;
       if strcmp(shape(i),'^')
          x_patch_tmp=[xcor(i) xcor(i)+s_patch/2 xcor(i)-s_patch/2 xcor(i)];
          y_patch_tmp=[ycor(i)+s_patch ycor(i)-s_patch/2 ycor(i)-s_patch/2 ycor(i)+s_patch];
          angle_rot=(heading(i)-90)*pi/180;
       elseif strcmp(shape(i),'d')
          x_patch_tmp=[xcor(i)+s_patch/2 xcor(i) xcor(i)-s_patch/2 xcor(i)];
          y_patch_tmp=[ycor(i) ycor(i)+1.2*s_patch ycor(i) ycor(i)-1.2*s_patch];
       elseif strcmp(shape(i),'s')
          x_patch_tmp=[xcor(i)-s_patch/2 xcor(i)+s_patch/2 xcor(i)+s_patch/2 xcor(i)-s_patch/2];
          y_patch_tmp=[ycor(i)+s_patch/2 ycor(i)+s_patch/2 ycor(i)-s_patch/2 ycor(i)-s_patch/2];
       end
       x_patch_tmp=x_patch_tmp-xcor(i);
       y_patch_tmp=y_patch_tmp-ycor(i);
       x_patch(i,:)=x_patch_tmp*cos(angle_rot)-y_patch_tmp*sin(angle_rot);
       y_patch(i,:)=y_patch_tmp*cos(angle_rot)+x_patch_tmp*sin(angle_rot);
       x_patch(i,:)=x_patch(i,:)+xcor(i);
       y_patch(i,:)=y_patch(i,:)+ycor(i);
    end

   return

end




